# import netflix.txt
netflix <- read.delim("~/Desktop/MA611/data/netflix.txt")
#---Creating a time-series object---#
sales=netflix[,3]
data=ts(sales,start=c(2000,1),frequency = 4)

#---Plotting it----#
autoplot(data) + ggtitle("Netflix Quartly Sales from 2000 to 2009")

#---Checking seasonality through polar-plot---#
ggseasonplot(data, year.labels=FALSE, continuous=TRUE, polar = TRUE)+
  ggtitle("Seasonality through polarmap for Netflix sales figures")

#---Checking seasonality through heatmap---#
Time.Stamp=seq(1,nrow(netflix),1)
data.modified=cbind(Time.Stamp,netflix)
ggplot(data.modified,aes(x = Time.Stamp, y = 1)) +
  geom_tile(aes(fill = sales)) +
  scale_fill_gradient2(low = "navy", mid = "yellow",
                       high = "red", midpoint=28) + ggtitle("Seasonality through heatmap for Netflix sales")+
  ylab("") + scale_y_discrete(expand=c(0,0))

# Classic Decomposition
#---Decomposing it and plotting the decomposition----#
dec=decompose(data)
plot(decompose(data))

#--Checking the strength of trend and seasonality for Netflix data---#
# FT
ft <- 1-var(dec$random,na.rm=TRUE)/var((dec$trend+dec$random),na.rm=TRUE)
# FS
fs <- 1-var(dec$random,na.rm=TRUE)/var((dec$seasonal+dec$random),na.rm=TRUE)
data.frame(cbind(ft,fs))
##          ft        fs
## 1 0.9990173 0.3788519
#---Now doing an stl decomposition---#
decomp=stl(data,t.window = 5,s.window="periodic", robust=TRUE)
decomp
##  Call:
##  stl(x = data, s.window = "periodic", t.window = 5, robust = TRUE)
## 
## Components
##           seasonal      trend     remainder
## 2000 Q1  1.3157587   3.854241  0.000000e+00
## 2000 Q2  0.9304170   7.362308 -1.142725e+00
## 2000 Q3 -0.6903739  10.870374  0.000000e+00
## 2000 Q4 -1.5558018  13.307308  1.638494e+00
## 2001 Q1  1.3157587  15.744241  0.000000e+00
## 2001 Q2  0.9304170  17.657308 -2.277245e-01
## 2001 Q3 -0.6903739  19.570374  0.000000e+00
## 2001 Q4 -1.5558018  24.392308 -1.216506e+00
## 2002 Q1  1.3157587  29.214241  0.000000e+00
## 2002 Q2  0.9304170  35.317308  1.122755e-01
## 2002 Q3 -0.6903739  41.420374  0.000000e+00
## 2002 Q4 -1.5558018  47.887308 -1.141506e+00
## 2003 Q1  1.3157587  54.354241  0.000000e+00
## 2003 Q2  0.9304170  63.622308 -1.362725e+00
## 2003 Q3 -0.6903739  72.890374  0.000000e+00
## 2003 Q4 -1.5558018  85.697308 -2.951506e+00
## 2004 Q1  1.3157587  98.504241  3.126388e-13
## 2004 Q2  0.9304170 118.779583  1.421085e-13
## 2004 Q3 -0.6903739 130.497692  1.060268e+01
## 2004 Q4 -1.5558018 142.215802  0.000000e+00
## 2005 Q1  1.3157587 152.657692 -1.523451e+00
## 2005 Q2  0.9304170 163.099583  0.000000e+00
## 2005 Q3 -0.6903739 178.827692 -5.397319e+00
## 2005 Q4 -1.5558018 194.555802  0.000000e+00
## 2006 Q1  1.3157587 216.487692  6.326549e+00
## 2006 Q2  0.9304170 238.419583  0.000000e+00
## 2006 Q3 -0.6903739 258.602692 -1.962319e+00
## 2006 Q4 -1.5558018 278.785802  0.000000e+00
## 2007 Q1  1.3157587 290.772692  1.323155e+01
## 2007 Q2  0.9304170 302.759583  0.000000e+00
## 2007 Q3 -0.6903739 303.337692 -8.677319e+00
## 2007 Q4 -1.5558018 303.915802  0.000000e+00
## 2008 Q1  1.3157587 320.297692  4.566549e+00
## 2008 Q2  0.9304170 336.679583  0.000000e+00
## 2008 Q3 -0.6903739 349.087692 -7.127319e+00
## 2008 Q4 -1.5558018 361.495802  0.000000e+00
## 2009 Q1  1.3157587 384.577692  8.206549e+00
## 2009 Q2  0.9304170 407.659583  0.000000e+00
## 2009 Q3 -0.6903739 426.702692 -2.892319e+00
## 2009 Q4 -1.5558018 445.745802  1.477929e-12
#---Plotting the stl decomposition---#
autoplot(decomp)+
  ggtitle("stl decomposition of Netflix data")

#---Now creating a lag-plot---#
gglagplot(data,lags=25,set.lags = 1:25)+
  ggtitle("Lag plots, Netflix Sales")

#----Now creating an ACF plot---#
ggAcf(data,lag.max = 30)+ggtitle("ACF plot for Netflix sales data")

#---Now creating a tapered ACF plot---#
ggtaperedacf(data,lag.max=25,calc.ci = T,level=95,nsim = 100)+
  ggtitle("Tapered ACF plot for Netflix sales data")

##################

#----Nomality checking-----#
##################
###################
###################
dec$random
##            Qtr1       Qtr2       Qtr3       Qtr4
## 2000         NA         NA  2.1181597  3.6124653
## 2001 -2.8598958 -1.4719792  0.6131597  0.5912153
## 2002 -3.1336458 -0.7769792  1.7819097  0.9174653
## 2003 -3.5436458 -2.2532292  1.0156597 -1.4062847
## 2004 -6.0211458  0.1130208 11.0781597  0.3812153
## 2005 -5.0636458 -1.8632292 -4.3780903 -1.3212847
## 2006  2.3388542 -1.1669792 -0.9668403  3.2937153
## 2007 10.9351042  3.6155208 -7.5755903 -4.8612847
## 2008  1.1526042  1.6767708 -6.0730903 -3.5937847
## 2009  3.8088542 -0.2594792         NA         NA
# overestimating
gghistogram(dec$random)
## Warning: Removed 4 rows containing non-finite values (stat_bin).

ggAcf(dec$random,lag.max=20)

#---Naive modeling---#
fitted(naive(data))
##        Qtr1   Qtr2   Qtr3   Qtr4
## 2000     NA   5.17   7.15  10.18
## 2001  13.39  17.06  18.36  18.88
## 2002  21.62  30.53  36.36  40.73
## 2003  45.19  55.67  63.19  72.20
## 2004  81.19  99.82 119.71 140.41
## 2005 140.66 152.45 164.03 172.74
## 2006 193.00 224.13 239.35 255.95
## 2007 277.23 305.32 303.69 293.97
## 2008 302.36 326.18 337.61 341.27
## 2009 359.94 394.10 408.59 423.12
#compare#
autoplot(data)+autolayer(fitted(naive(data)))+ggtitle("Tracking real and naively predicted Netflix sales")
## Warning: Removed 1 rows containing missing values (geom_path).

#--extract residuals--#
naive.model.resid=residuals(naive(data))
gghistogram(naive.model.resid) ## underestimating
## Warning: Removed 1 rows containing non-finite values (stat_bin).

ggAcf(naive.model.resid)

#--another way of extracting residuals---#
diff(sales)
##  [1]  1.98  3.03  3.21  3.67  1.30  0.52  2.74  8.91  5.83  4.37  4.46 10.48
## [13]  7.52  9.01  8.99 18.63 19.89 20.70  0.25 11.79 11.58  8.71 20.26 31.13
## [25] 15.22 16.60 21.28 28.09 -1.63 -9.72  8.39 23.82 11.43  3.66 18.67 34.16
## [37] 14.49 14.53 21.07
#--checking normality---#
gghistogram(naive.model.resid)
## Warning: Removed 1 rows containing non-finite values (stat_bin).

ggAcf(naive.model.resid,lag.max=20)

autoplot(naive.model.resid)+ggtitle("Naive residuals showing reasonable stationarity")

netflix.ts <- data
#---State-space ETS modeling---#
########
intuitive.model=ets(netflix.ts,model="AAN")
intuitive.model
## ETS(A,A,N) 
## 
## Call:
##  ets(y = netflix.ts, model = "AAN") 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.1896 
## 
##   Initial states:
##     l = -0.1455 
##     b = 3.5045 
## 
##   sigma:  9.3247
## 
##      AIC     AICc      BIC 
## 331.9543 333.7190 340.3986
accuracy(intuitive.model)
##                    ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
## Training set 1.749013 8.846204 5.942018 1.525494 5.891768 0.1309047 0.2066377
#autoplot(intuitive.model)
autoplot(fitted(intuitive.model))+autolayer(netflix.ts)+ggtitle("Comparing the intuitive ETS(AAN) fitted values to the actuals, Netflix")+ylab("Observed/Fitted")

forecast(intuitive.model,h=10)
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2010 Q1       460.9564 449.0063 472.9065 442.6803 479.2325
## 2010 Q2       477.7233 459.1531 496.2935 449.3227 506.1239
## 2010 Q3       494.4902 469.6622 519.3182 456.5190 532.4614
## 2010 Q4       511.2571 480.1475 542.3667 463.6791 558.8352
## 2011 Q1       528.0240 490.4841 565.5640 470.6116 585.4364
## 2011 Q2       544.7909 500.6207 588.9612 477.2384 612.3435
## 2011 Q3       561.5579 510.5347 612.5811 483.5246 639.5911
## 2011 Q4       578.3248 520.2165 636.4331 489.4558 667.1937
## 2012 Q1       595.0917 529.6631 660.5203 495.0273 695.1561
## 2012 Q2       611.8586 538.8750 684.8422 500.2398 723.4774
autoplot(forecast(intuitive.model,h=10))+ggtitle("Extracting forecasts from the intuitive ETS(AAN), Netflix")+ylab("Values")

best.model=ets(netflix.ts)
best.model
## ETS(M,A,M) 
## 
## Call:
##  ets(y = netflix.ts) 
## 
##   Smoothing parameters:
##     alpha = 0.9997 
##     beta  = 0.84 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 2.2848 
##     b = 2.9875 
##     s = 0.9714 0.9843 1.0124 1.0319
## 
##   sigma:  0.062
## 
##      AIC     AICc      BIC 
## 301.2328 307.2328 316.4327
accuracy(best.model)
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.6882066 6.837355 4.806047 0.5518281 3.945826 0.1058789
##                    ACF1
## Training set -0.1473851
#autoplot(best.model)
autoplot(fitted(best.model))+autolayer(netflix.ts)+ggtitle("Comparing the best ETS(MAM) fitted values to the actuals, Netflix")+ylab("Observed/Fitted")

forecast(best.model,h=10)
##         Point Forecast     Lo 80     Hi 80      Lo 95     Hi 95
## 2010 Q1       499.8327 460.14173  539.5236  439.13062  560.5347
## 2010 Q2       517.7683 435.05494  600.4816  391.26913  644.2674
## 2010 Q3       530.0811 397.43001  662.7322  327.20875  732.9535
## 2010 Q4       549.4192 357.91313  740.9252  256.53595  842.3024
## 2011 Q1       611.6243 334.41903  888.8296  187.67541 1035.5732
## 2011 Q2       627.4393 273.94448  980.9341   86.81563 1168.0629
## 2011 Q3       636.7137 204.61719 1068.8101  -24.12087 1297.5482
## 2011 Q4       654.6498 131.70031 1177.5993 -145.13242 1454.4320
## 2012 Q1       723.4168  54.94848 1391.8851 -298.91727 1745.7509
## 2012 Q2       737.1111 -40.20095 1514.4232 -451.68510 1925.9073
autoplot(forecast(best.model,h=10))+ggtitle("Extracting forecasts from the best ETS(MAM), Netflix")+ylab("Values")

##########
##########
#---Comparing different models---#
##########
##########

illogical.model=ets(netflix.ts,model="ANN")
illogical.model
## ETS(A,N,N) 
## 
## Call:
##  ets(y = netflix.ts, model = "ANN") 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
## 
##   Initial states:
##     l = 5.2444 
## 
##   sigma:  14.9002
## 
##      AIC     AICc      BIC 
## 367.6136 368.2803 372.6802
accuracy(illogical.model)
##                    ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
## Training set 10.97469 14.52295 11.54577 10.11385 10.37791 0.2543573 0.4209748
autoplot(fitted(illogical.model))+autolayer(netflix.ts)+ggtitle("Comparing the illogical ETS(ANN) fitted values to the actuals, Netflix")+ylab("Observed/Fitted")

naive.model=naive(netflix.ts)
naive.model
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2010 Q1         444.19 425.3424 463.0376 415.3651 473.0149
## 2010 Q2         444.19 417.5354 470.8446 403.4254 484.9546
## 2010 Q3         444.19 411.5450 476.8350 394.2637 494.1163
## 2010 Q4         444.19 406.4948 481.8852 386.5401 501.8399
## 2011 Q1         444.19 402.0454 486.3346 379.7355 508.6445
## 2011 Q2         444.19 398.0230 490.3570 373.5836 514.7964
## 2011 Q3         444.19 394.3239 494.0561 367.9264 520.4536
## 2011 Q4         444.19 390.8809 497.4991 362.6607 525.7193
## 2012 Q1         444.19 387.6471 500.7329 357.7152 530.6648
## 2012 Q2         444.19 384.5886 503.7914 353.0375 535.3425
accuracy(naive.model)
##                    ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
## Training set 11.25692 14.70688 11.83897 10.40919 10.60627 0.2608166 0.4067258
autoplot(fitted(naive.model))+autolayer(netflix.ts)+ggtitle("Comparing the naive fitted values to the actuals, Netflix")+ylab("Observed/Fitted")

forecast(naive.model,h=10)
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2010 Q1         444.19 425.3424 463.0376 415.3651 473.0149
## 2010 Q2         444.19 417.5354 470.8446 403.4254 484.9546
## 2010 Q3         444.19 411.5450 476.8350 394.2637 494.1163
## 2010 Q4         444.19 406.4948 481.8852 386.5401 501.8399
## 2011 Q1         444.19 402.0454 486.3346 379.7355 508.6445
## 2011 Q2         444.19 398.0230 490.3570 373.5836 514.7964
## 2011 Q3         444.19 394.3239 494.0561 367.9264 520.4536
## 2011 Q4         444.19 390.8809 497.4991 362.6607 525.7193
## 2012 Q1         444.19 387.6471 500.7329 357.7152 530.6648
## 2012 Q2         444.19 384.5886 503.7914 353.0375 535.3425
autoplot(forecast(naive.model,h=10))+ggtitle("Extracting forecasts from the naive model, Netflix")+ylab("Values")

stl.randomwalk.model=stl(netflix.ts,t.window = 5,s.window="periodic", robust=TRUE)
stl.randomwalk.model
##  Call:
##  stl(x = netflix.ts, s.window = "periodic", t.window = 5, robust = TRUE)
## 
## Components
##           seasonal      trend     remainder
## 2000 Q1  1.3157587   3.854241  0.000000e+00
## 2000 Q2  0.9304170   7.362308 -1.142725e+00
## 2000 Q3 -0.6903739  10.870374  0.000000e+00
## 2000 Q4 -1.5558018  13.307308  1.638494e+00
## 2001 Q1  1.3157587  15.744241  0.000000e+00
## 2001 Q2  0.9304170  17.657308 -2.277245e-01
## 2001 Q3 -0.6903739  19.570374  0.000000e+00
## 2001 Q4 -1.5558018  24.392308 -1.216506e+00
## 2002 Q1  1.3157587  29.214241  0.000000e+00
## 2002 Q2  0.9304170  35.317308  1.122755e-01
## 2002 Q3 -0.6903739  41.420374  0.000000e+00
## 2002 Q4 -1.5558018  47.887308 -1.141506e+00
## 2003 Q1  1.3157587  54.354241  0.000000e+00
## 2003 Q2  0.9304170  63.622308 -1.362725e+00
## 2003 Q3 -0.6903739  72.890374  0.000000e+00
## 2003 Q4 -1.5558018  85.697308 -2.951506e+00
## 2004 Q1  1.3157587  98.504241  3.126388e-13
## 2004 Q2  0.9304170 118.779583  1.421085e-13
## 2004 Q3 -0.6903739 130.497692  1.060268e+01
## 2004 Q4 -1.5558018 142.215802  0.000000e+00
## 2005 Q1  1.3157587 152.657692 -1.523451e+00
## 2005 Q2  0.9304170 163.099583  0.000000e+00
## 2005 Q3 -0.6903739 178.827692 -5.397319e+00
## 2005 Q4 -1.5558018 194.555802  0.000000e+00
## 2006 Q1  1.3157587 216.487692  6.326549e+00
## 2006 Q2  0.9304170 238.419583  0.000000e+00
## 2006 Q3 -0.6903739 258.602692 -1.962319e+00
## 2006 Q4 -1.5558018 278.785802  0.000000e+00
## 2007 Q1  1.3157587 290.772692  1.323155e+01
## 2007 Q2  0.9304170 302.759583  0.000000e+00
## 2007 Q3 -0.6903739 303.337692 -8.677319e+00
## 2007 Q4 -1.5558018 303.915802  0.000000e+00
## 2008 Q1  1.3157587 320.297692  4.566549e+00
## 2008 Q2  0.9304170 336.679583  0.000000e+00
## 2008 Q3 -0.6903739 349.087692 -7.127319e+00
## 2008 Q4 -1.5558018 361.495802  0.000000e+00
## 2009 Q1  1.3157587 384.577692  8.206549e+00
## 2009 Q2  0.9304170 407.659583  0.000000e+00
## 2009 Q3 -0.6903739 426.702692 -2.892319e+00
## 2009 Q4 -1.5558018 445.745802  1.477929e-12
stl.forecast=forecast(stl.randomwalk.model, h=10,method="naive",robust=FALSE)
stl.randomwalk.model$time.series[,1]+stl.randomwalk.model$time.series[,2]
##            Qtr1       Qtr2       Qtr3       Qtr4
## 2000   5.170000   8.292725  10.180000  11.751506
## 2001  17.060000  18.587725  18.880000  22.836506
## 2002  30.530000  36.247725  40.730000  46.331506
## 2003  55.670000  64.552725  72.200000  84.141506
## 2004  99.820000 119.710000 129.807319 140.660000
## 2005 153.973451 164.030000 178.137319 193.000000
## 2006 217.803451 239.350000 257.912319 277.230000
## 2007 292.088451 303.690000 302.647319 302.360000
## 2008 321.613451 337.610000 348.397319 359.940000
## 2009 385.893451 408.590000 426.012319 444.190000
autoplot(stl.randomwalk.model$time.series[,1]+stl.randomwalk.model$time.series[,2])+autolayer(netflix.ts)+ggtitle("Comparing the stl signal values to the actuals, Netflix")+ylab("Observed/stl.Signal")

autoplot(stl.forecast)+ggtitle("Extracting forecasts from the stl model, Netflix")+ylab("Values")

accuracy(stl.forecast)
##                    ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
## Training set 11.33055 14.36951 11.80973 10.93826 11.10056 0.2601723 0.4861571
###################
#---Checking residual assumptions------#
###################

checkresiduals(forecast(naive.model))

## 
##  Ljung-Box test
## 
## data:  Residuals from Naive method
## Q* = 19.518, df = 8, p-value = 0.01232
## 
## Model df: 0.   Total lags used: 8
checkresiduals(forecast(illogical.model))

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,N,N)
## Q* = 20.246, df = 6, p-value = 0.002504
## 
## Model df: 2.   Total lags used: 8
checkresiduals(forecast(intuitive.model))

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,A,N)
## Q* = 28.556, df = 4, p-value = 9.62e-06
## 
## Model df: 4.   Total lags used: 8
checkresiduals(forecast(best.model))

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(M,A,M)
## Q* = 31.454, df = 3, p-value = 6.82e-07
## 
## Model df: 8.   Total lags used: 11
###############
#--Forced damping--#
###############

damped.model=ets(netflix.ts,model="AAN",damped = TRUE)
damped.model
## ETS(A,Ad,N) 
## 
## Call:
##  ets(y = netflix.ts, model = "AAN", damped = TRUE) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
##     beta  = 0.2737 
##     phi   = 0.9799 
## 
##   Initial states:
##     l = 4.6853 
##     b = 1.5605 
## 
##   sigma:  9.5868
## 
##      AIC     AICc      BIC 
## 335.0446 337.5900 345.1779
accuracy(damped.model)
##                    ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
## Training set 2.074655 8.967594 6.024618 2.425019 5.713945 0.1327244 0.1601022
forecast(damped.model,h=100)
##         Point Forecast     Lo 80     Hi 80       Lo 95     Hi 95
## 2010 Q1       460.8196 448.53371  473.1056  442.029935  479.6093
## 2010 Q2       477.1161 457.27494  496.9572  446.771673  507.4605
## 2010 Q3       493.0855 465.74624  520.4247  451.273730  534.8972
## 2010 Q4       508.7344 473.65972  543.8090  455.092327  562.3764
## 2011 Q1       524.0692 480.95426  567.1841  458.130605  590.0077
## 2011 Q2       539.0962 487.62629  590.5662  460.379751  617.8127
## 2011 Q3       553.8217 493.69122  613.9522  461.860065  645.7833
## 2011 Q4       568.2516 499.17112  637.3322  462.602104  673.9012
## 2012 Q1       582.3920 504.09007  660.6939  462.639547  702.1444
## 2012 Q2       596.2485 508.47222  684.0248  462.006260  730.4908
## 2012 Q3       609.8270 512.34100  707.3129  460.735048  758.9189
## 2012 Q4       623.1329 515.71877  730.5470  458.857159  787.4087
## 2013 Q1       636.1718 518.62672  753.7169  456.402115  815.9415
## 2013 Q2       648.9490 521.08489  776.8131  453.397720  844.5003
## 2013 Q3       661.4698 523.11220  799.8274  449.870121  873.0694
## 2013 Q4       673.7393 524.72651  822.7520  445.843913  901.6346
## 2014 Q1       685.7625 525.94469  845.5803  441.342247  930.1828
## 2014 Q2       697.5445 526.78274  868.3062  436.386941  958.7020
## 2014 Q3       709.0899 527.25578  890.9241  430.998580  987.1813
## 2014 Q4       720.4037 527.37817  913.4293  425.196618 1015.6108
## 2015 Q1       731.4904 527.16357  935.8173  418.999459 1043.9814
## 2015 Q2       742.3546 526.62494  958.0843  412.424541 1072.2847
## 2015 Q3       753.0008 525.77466  980.2269  405.488407 1100.5132
## 2015 Q4       763.4333 524.62451 1002.2421  398.206769 1128.6598
## 2016 Q1       773.6564 523.18574 1024.1271  390.594569 1156.7183
## 2016 Q2       783.6744 521.46912 1045.8797  382.666031 1184.6827
## 2016 Q3       793.4913 519.48492 1067.4977  374.434708 1212.5479
## 2016 Q4       803.1112 517.24299 1088.9793  365.913528 1240.3088
## 2017 Q1       812.5380 514.75279 1110.3232  357.114833 1267.9611
## 2017 Q2       821.7756 512.02335 1131.5278  348.050416 1295.5008
## 2017 Q3       830.8278 509.06336 1152.5923  338.731553 1322.9241
## 2017 Q4       839.6984 505.88118 1173.5156  329.169033 1350.2277
## 2018 Q1       848.3909 502.48481 1194.2970  319.373189 1377.4086
## 2018 Q2       856.9090 498.88196 1214.9359  309.353921 1404.4640
## 2018 Q3       865.2561 495.08006 1235.4321  299.120719 1431.3914
## 2018 Q4       873.4357 491.08624 1255.7851  288.682687 1458.1886
## 2019 Q1       881.4511 486.90739 1275.9948  278.048565 1484.8536
## 2019 Q2       889.3057 482.55012 1296.0612  267.226742 1511.3846
## 2019 Q3       897.0026 478.02082 1315.9843  256.225278 1537.7799
## 2019 Q4       904.5450 473.32566 1335.7644  245.051920 1564.0381
## 2020 Q1       911.9361 468.47059 1355.4016  233.714116 1590.1581
## 2020 Q2       919.1789 463.46132 1374.8964  222.219028 1616.1387
## 2020 Q3       926.2762 458.30340 1394.2491  210.573549 1641.9789
## 2020 Q4       933.2312 453.00219 1413.4602  198.784312 1667.6781
## 2021 Q1       940.0465 447.56283 1432.5303  186.857704 1693.2354
## 2021 Q2       946.7251 441.99033 1451.4599  174.799873 1718.6504
## 2021 Q3       953.2697 436.28951 1470.2498  162.616745 1743.9226
## 2021 Q4       959.6829 430.46503 1488.9007  150.314027 1769.0517
## 2022 Q1       965.9674 424.52140 1507.4133  137.897222 1794.0375
## 2022 Q2       972.1257 418.46298 1525.7884  125.371631 1818.8798
## 2022 Q3       978.1605 412.29399 1544.0270  112.742369 1843.5786
## 2022 Q4       984.0741 406.01852 1562.1297  100.014366 1868.1339
## 2023 Q1       989.8691 399.64051 1580.0977   87.192378 1892.5458
## 2023 Q2       995.5478 393.16380 1597.9317   74.280993 1916.8145
## 2023 Q3      1001.1125 386.59208 1615.6328   61.284638 1940.9403
## 2023 Q4      1006.5655 379.92893 1633.2020   48.207584 1964.9234
## 2024 Q1      1011.9090 373.17783 1650.6403   35.053954 1988.7641
## 2024 Q2      1017.1454 366.34215 1667.9486   21.827728 2012.4630
## 2024 Q3      1022.2766 359.42513 1685.1281    8.532746 2036.0205
## 2024 Q4      1027.3049 352.42993 1702.1798   -4.827282 2059.4370
## 2025 Q1      1032.2322 345.35962 1719.1048  -18.248777 2082.7132
## 2025 Q2      1037.0607 338.21714 1735.9042  -31.728280 2105.8496
## 2025 Q3      1041.7922 331.00538 1752.5791  -45.262452 2128.8469
## 2025 Q4      1046.4288 323.72711 1769.1305  -58.848070 2151.7057
## 2026 Q1      1050.9724 316.38502 1785.5597  -72.482019 2174.4267
## 2026 Q2      1055.4247 308.98174 1801.8677  -86.161290 2197.0107
## 2026 Q3      1059.7877 301.51980 1818.0556  -99.882979 2219.4584
## 2026 Q4      1064.0631 294.00164 1834.1246 -113.644279 2241.7705
## 2027 Q1      1068.2527 286.42966 1850.0758 -127.442479 2263.9480
## 2027 Q2      1072.3583 278.80616 1865.9104 -141.274960 2285.9915
## 2027 Q3      1076.3814 271.13338 1881.6295 -155.139194 2307.9021
## 2027 Q4      1080.3238 263.41348 1897.2342 -169.032736 2329.6804
## 2028 Q1      1084.1871 255.64858 1912.7256 -182.953225 2351.3274
## 2028 Q2      1087.9729 247.84071 1928.1050 -196.898381 2372.8441
## 2028 Q3      1091.6826 239.99185 1943.3734 -210.866002 2394.2312
## 2028 Q4      1095.3179 232.10393 1958.5319 -224.853960 2415.4898
## 2029 Q1      1098.8803 224.17880 1973.5818 -238.860199 2436.6208
## 2029 Q2      1102.3711 216.21826 1988.5240 -252.882734 2457.6250
## 2029 Q3      1105.7919 208.22408 2003.3598 -266.919648 2478.5035
## 2029 Q4      1109.1441 200.19794 2018.0902 -280.969087 2499.2573
## 2030 Q1      1112.4290 192.14149 2032.7164 -295.029264 2519.8872
## 2030 Q2      1115.6479 184.05633 2047.2395 -309.098452 2540.3943
## 2030 Q3      1118.8022 175.94401 2061.6605 -323.174982 2560.7795
## 2030 Q4      1121.8933 167.80603 2075.9805 -337.257243 2581.0438
## 2031 Q1      1124.9223 159.64385 2090.2007 -351.343681 2601.1882
## 2031 Q2      1127.8905 151.45888 2104.3221 -365.432793 2621.2138
## 2031 Q3      1130.7991 143.25248 2118.3458 -379.523132 2641.1214
## 2031 Q4      1133.6494 135.02599 2132.2728 -393.613298 2660.9121
## 2032 Q1      1136.4424 126.78070 2146.1042 -407.701940 2680.5868
## 2032 Q2      1139.1794 118.51785 2159.8410 -421.787756 2700.1466
## 2032 Q3      1141.8615 110.23866 2173.4844 -435.869490 2719.5925
## 2032 Q4      1144.4898 101.94430 2187.0352 -449.945927 2738.9254
## 2033 Q1      1147.0653  93.63591 2200.4946 -464.015898 2758.1464
## 2033 Q2      1149.5891  85.31460 2213.8635 -478.078275 2777.2564
## 2033 Q3      1152.0622  76.98143 2227.1430 -492.131970 2796.2564
## 2033 Q4      1154.4857  68.63744 2240.3341 -506.175933 2815.1474
## 2034 Q1      1156.8606  60.28364 2253.4376 -520.209153 2833.9304
## 2034 Q2      1159.1878  51.92100 2266.4547 -534.230657 2852.6064
## 2034 Q3      1161.4684  43.55048 2279.3862 -548.239503 2871.1762
## 2034 Q4      1163.7031  35.17298 2292.2332 -562.234788 2889.6410
autoplot(forecast(damped.model,h=100))+ggtitle("Extracting forecasts from trend-damped AAN, Netflix")

checkresiduals(damped.model)

## 
##  Ljung-Box test
## 
## data:  Residuals from ETS(A,Ad,N)
## Q* = 31.608, df = 3, p-value = 6.33e-07
## 
## Model df: 5.   Total lags used: 8
##############
##############
##############
################
#---On simulations---#
################
simulate(intuitive.model)
##           Qtr1      Qtr2      Qtr3      Qtr4
## 2010  483.7831  513.3689  528.2873  553.9979
## 2011  575.1888  592.3600  608.9938  646.9377
## 2012  668.4516  707.5872  741.1059  769.0410
## 2013  800.5554  838.7958  880.4491  918.6317
## 2014  950.1046  969.6539 1008.1170 1038.3941
## 2015 1069.0269 1108.6436 1117.9986 1150.9169
## 2016 1191.8962 1233.6088 1274.5044 1314.2884
## 2017 1365.2982 1409.7686 1454.3235 1493.7692
## 2018 1532.4007 1582.7033 1666.0219 1711.8642
## 2019 1758.1163 1805.5605 1855.5084 1903.8551
autoplot(netflix.ts)+
  autolayer(simulate(illogical.model,future=F),series = "simANN1")+
  autolayer(simulate(illogical.model,future=F),series = "simANN2")+
  autolayer(simulate(damped.model,future=F),series = "simdamped1")+
  autolayer(simulate(damped.model,future=F),series = "simdamped2")+ylab("Real and simulated values")+
  ggtitle("Simulating the past, Netflix")

autoplot(netflix.ts)+
  autolayer(simulate(illogical.model,future=T),series = "simANN1")+
  autolayer(simulate(illogical.model,future=T),series = "simANN2")+
  autolayer(simulate(damped.model,future=T),series = "simdamped1")+
  autolayer(simulate(damped.model,future=T),series = "simdamped2")+ylab("Real and simulated values")+
  ggtitle("Simulating the future, Netflix")

######
#####
MAM1=simulate(best.model,future=F)
MAM2=simulate(best.model,future=F)
MAM3=simulate(best.model,future=F)
MAM4=simulate(best.model,future=F)

AAN1=simulate(intuitive.model,future=F)
AAN2=simulate(intuitive.model,future=F)
AAN3=simulate(intuitive.model,future=F)
AAN4=simulate(intuitive.model,future=F)

ANN1=simulate(illogical.model,future=F)
ANN2=simulate(illogical.model,future=F)
ANN3=simulate(illogical.model,future=F)
ANN4=simulate(illogical.model,future=F)

#Naive1=simulate(naive.model,future=F)
#Naive2=simulate(naive.model,future=F)
#Naive3=simulate(naive.model,future=F)
#Naive4=simulate(naive.model,future=F)

data=ts(cbind(MAM1,MAM2,MAM3,MAM4,AAN1,AAN2,AAN3,AAN4,ANN1,ANN2,ANN3,ANN4,netflix.ts),start=c(2000,1),frequency = 4)


dissimilarity=diss(data,METHOD="COR")
dissimilarity
##                 MAM1      MAM2      MAM3      MAM4      AAN1      AAN2
## MAM2       0.6151954                                                  
## MAM3       1.9673322 1.8937494                                        
## MAM4       0.6765662 1.1148670 1.8872012                              
## AAN1       1.9188661 1.8390623 0.4759949 1.8092863                    
## AAN2       1.8748312 1.9462609 0.6311500 1.6014022 0.6472225          
## AAN3       1.8994251 1.7473616 0.5474444 1.8801718 0.3000228 0.8747049
## AAN4       1.9528903 1.9165466 0.4119907 1.7671558 0.3131599 0.4190454
## ANN1       1.5827130 1.2199543 1.2788737 1.7413233 1.1375033 1.5457544
## ANN2       0.9274408 1.0357490 1.8183985 1.1382533 1.8355100 1.8030943
## ANN3       1.1177687 1.4164630 1.6928040 0.9588245 1.7725410 1.4929428
## ANN4       0.5351675 0.7203757 1.9533664 0.8419682 1.9920761 1.9051435
## netflix.ts 1.9473603 1.9515145 0.4231992 1.7337572 0.4296210 0.3115663
##                 AAN3      AAN4      ANN1      ANN2      ANN3      ANN4
## MAM2                                                                  
## MAM3                                                                  
## MAM4                                                                  
## AAN1                                                                  
## AAN2                                                                  
## AAN3                                                                  
## AAN4       0.5481173                                                  
## ANN1       0.9916845 1.2972458                                        
## ANN2       1.8125665 1.8340317 1.5153817                              
## ANN3       1.8210869 1.6702093 1.8602124 1.0663360                    
## ANN4       1.9667951 1.9765064 1.5989999 0.8151843 1.0100359          
## netflix.ts 0.6609106 0.1771690 1.3879294 1.8288415 1.6104559 1.9613553
hc.dpred <- hclust(dissimilarity)
plot(hc.dpred,main="Cluster dendogram, Netflix, Correlation distance")